Accommodating detection limits of multiple exposures in environmental mixture analyses: an overview of statistical approaches

Background Identifying the impact of environmental mixtures on human health is an important topic. However, such studies face challenges when exposure measurements lie below limit of detection (LOD). While various approaches for accommodating a single exposure subject to LOD have been used, their impact on mixture analysis has not been thoroughly investigated. Our study aims to understand the impact of five popular LOD accommodation approaches on mixture analysis results with multiple exposures subject to LOD, including omitting subjects with any exposures below LOD (complete case analysis); single imputations by LOD/\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{2}$$\end{document}2, and by estimates from a censored accelerated failure time (AFT) model; and multiple imputation (MI) with or without truncation based on LOD. Methods In extensive simulation studies with high-dimensional and highly correlated exposures and a continuous health outcome, we examined the performance of each LOD approach on three mixture analysis methods: elastic net regression, weighted quantile sum regression (WQS) and Bayesian kernel machine regression (BKMR). We further analyzed data from the National Health and Nutrition Examination Survey (NHANES) on how persistent organic pollutants (POPs) influenced leukocyte telomere length (LTL). Results Complete case analysis was inefficient and could result in severe bias for some mixture methods. Imputation by LOD/\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{2}$$\end{document}2 showed unstable performance across mixture methods. Conventional MI was associated with consistent mild biases, which can be reduced by using a truncated distribution for imputation. Estimating censored values by AFT models had a minimal impact on the results. In the NHANES analysis, imputation by LOD/\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{2}$$\end{document}2, truncated MI and censored AFT models performed similarly, with a positive overall effect of POPs on LTL while PCB126, PCB169 and furan 2,3,4,7,8-pncdf being the most important exposures. Conclusions Our study favored using truncated MI and censored AFT models to accommodate values below LOD for the stability of downstream mixture analysis. Supplementary Information The online version contains supplementary material available at 10.1186/s12940-024-01088-w.


Background
Environmental exposures to chemical, biological, or physical substances found in air, water, food, or soil are common during the human life course [1][2][3].These highdimensional and highly correlated exposures can act synergistically or antagonistically on human health [4,5].Studying individual exposures only addresses their marginal effects, without accounting for others, which can result in misleading conclusions about effects of the whole mixture [6,7].
Several popular modeling approaches exist to analyze complex environmental mixtures, including but not limited to regularized regressions, weighted quantile sum regression (WQS) [8,9], and Bayesian kernel machine regression (BKMR) [10,11].Briefly, regularized regressions such as elastic net regression [12] and lasso (least absolute shrinkage and selection operator) [13] can be used to identify the relative importance of driver(s) in the mixture through variable selection [14][15][16].WQS derives a one-dimensional weighted sum score of the exposures with a linear relationship to a continuous health outcome under the assumption that all exposure effects are in the same direction.WQS has been generalized to several types of outcomes [17] and is widely used in practice [16,[18][19][20].BKMR is a Bayesian nonparametric method to handle complex nonlinear relationships between exposure mixtures and continuous, binary, and time-to-event outcomes [10,21].It has been widely used in mixtures studies, due to its flexibility and abundant visualization tools [16,22].Details of these methods are illustrated in Appendix A.
Environmental health studies often face challenges of exposure values below limit of detection (LOD) (i.e., leftcensored).All the above-mentioned mixture methods assume accurate measurements of exposures, thus some procedure for accommodating LOD is needed before applying these mixture methods.For example, with data from the National Health and Nutrition Examination Survey (NHANES) 2001-2002 cycle, Gibson et al. [23] investigated the relationship between persistent organic pollutants (POPs) and leukocyte telomere length (LTL), a biomarker associated with chronic diseases [24][25][26][27] and dioxin-associated cancers [28][29][30][31].Among the 34 POPs with 1.4% to 99.9% of values below LOD [32], Gibson et al. [23] restricted their analysis to the 18 POPs with less than 40% of values below LOD and imputed all values below LOD by LOD/ √ 2 .However, it is unclear how this imputation influenced the analysis results.
Recovering the true effects of environmental mixtures, where multiple exposures are subject to different proportions of values below LOD, is thus an important problem to address.Several approaches for accommodating values below LOD for a single exposure have been used in practice, including complete case analysis by omitting subjects with any measured values below LOD, single imputation by LOD/ √ 2 , and multiple imputa- tion.Ortega-Villa et al. [33] empirically compared these approaches in an environmental study setting with a binary outcome and a single exposure.However, the impact of these LOD approaches on downstream mixture analysis results has not been thoroughly investigated in settings where multiple exposures within the highdimensional and highly correlated exposure mixtures are subject to LOD.
In this manuscript, we aim to understand the impact of five popular approaches for accommodating LOD: complete case analysis; single imputation of values below LOD by LOD/ √ 2 and by estimates from censored accel- erated failure time (AFT) models; and multiple imputation (MI) with and without LOD-based truncation.We conducted extensive simulation studies to examine their influences on three popular mixture methods, including elastic net regression, WQS and BKMR, after applying the above-mentioned LOD approaches.We also re-analyzed the 2001-2002 NHANES dataset as described in Gibson et al. [23], to illustrate how different ways of handling LOD can impact the identification of associations between the POPs and LTL.Through these simulated and real data examples, we would like to draw readers' attention to carefully choose LOD accommodation approaches for mixture analysis, rather than recommending one approach as the gold standard.

LOD accommodation approaches
Here we give a brief review of the five LOD accommodation approaches.Complete case analysis only includes subjects whose exposure values are all above LOD.In theory, this approach provides unbiased results for linear regressions when the missingness only depends on the exposures [34,35].However, its performance may be unstable in practice with reduced sample sizes [36,37].An alternative approach is to replace values below LOD with a pre-specified constant value such as LOD, LOD/2, or LOD/ √ 2 based on the observed exposure distribution [37][38][39].In this study we chose LOD/ √ 2 , which is widely used for log-normally distributed (or right skewed) chemical exposures.This approach is popular due to its simplicity, but results may be biased when the distribution of values below LOD is not centered on the substitution value [37,38,40].
Chen et al. [41] recently proposed a new approach using multivariate accelerated failure time (AFT) regressions to model multiple left-censored chemicals through baseline covariates, which is a flexible approach specialized to handle censored outcomes with mild assumptions about the joint distribution.Note that this is an extension of the approach proposed in Kong and Nan [42] from a single exposure subject to LOD to multiple exposures.Due to simultaneously fitting multiple AFT models, it allows one to specifically account for the correlations between chemicals through shared baseline covariates and correlation between error terms.The originally proposed approach in Chen et al. [41] allows one to simultanueous model the exposures and health outcomes with efficiency gain.However, due to practical considerations, we only adopted the first part of this approach with the multivariate AFT model to conduct a single imputation for simplicity.The details of this approach are described in Appendix B.
Lastely, instead of single imputation approaches described above, multiple imputation (MI) also has been widely used by treating values below LOD as missing, and then imputing with models such as Bayesian linear regression or linear regression with bootstrap samples [35,43].MI generates multiple datasets (e.g., 5 or 10) for downstream analysis and combine analysis results using the Rubin's rule [44].In this study we chose to use the bootstrap linear regression implemented in the R 'mice' package [45] due to its superior performance in the settings we investigated.However, conventional MI does not guarantee that imputed values are below LOD.Thus, we improved it by truncating the estimated normal distribution at LOD to ensure all imputed values are in the correct range, and named this approach as truncated MI.The details of conventional and truncated MI are described in Appendix C.

Simulation settings
We conducted extensive simulations to empirically evaluate the impact of LOD accommodation approaches on three popular downstream mixture analysis methods, including elastic net regression, WQS and BKMR, under various settings.First, covariates X = (1, X 1 , X 2 ) T were independently generated from X 1 ∼ Bern(p = 0.5) and X 2 ∼ N (1, 1) .Given that environmental exposures are commonly highly correlated, right-skewed and associated through covariates, a mixture of p = 10 exposures Z = (Z 1 , . . ., Z 10 ) T was generated from a multivariate linear regression model with covariates X and log link, that is, , where R 1 and R 2 are 3 × 3 correlation matrices with all off-diagonal entries as 0.25 and 0.75, respectively, and R 3 is a 4 × 4 correlation matrix with all off-diagonal entries as 0.5.Through this formulation, we imposed correlations between exposures through two sources: shared covariate effects X , where the correlations are governed by η , and correlation between error terms through off-diagonal entries in .By the group structure in (i.e., {Z 1 , Z 2 , Z 3 } for group 1, {Z 4 , Z 5 , Z 6 } for group 2, and {Z 7 , Z 8 , Z 9 , Z 10 } for group 3), we allowed a higher within-group correlation than between-group correlations.We varied σ = 1/2 and 1/8 for moderate and high correlations within the groups, respectively (see Figure S1 for Spearman correlation coefficients between simulated variables).Because Z were right-skewed, we generated outcome Y under a linear regression with the log-transformed Z log , as in many envi- ronmental health studies, that is, Y = β T Z log + α T X + ǫ, where ǫ ∼ N (0, 2) .With a sample size of 500, we fixed α = (1, 1, 1) , and varied β and percent of value below LOD in various scenarios as follows.
-Scenario 1.We set β = (1.0,0.8, 0.0, 0.6, 0.4, 0.0, 0.2, 0.1, 0.0, 0.0) T to reflect the relative importance of these exposures, and assumed that Z 2 , Z 3 , Z 5 , Z 7 , and Z 9 have approxi- mately 30% of values below LOD.-Scenario 2. Z 2 is assumed to have approximately 70% of values below LOD, while all the other settings are the same as in Scenario 1.In this scenario, we handled Z 2 in two ways that are widely used in practice: (i) Z 2 was completely excluded from the analysis (Scenario 2A), and (ii) an indicator variable of whether Z 2 is above the LOD was used (Scenario 2B), while the other exposures subject to LOD were handled with the above-mentioned approaches.This scenario allows us to understand how to handle an exposure with a high percent of values below LOD.-Scenario 3. We generated all the exposures as in Scenario 1, but we re-generated a new Z 2 from Unif (0, LOD) if the original Z 2 was below LOD.This essentially resulted in the marginal distribution of Z 2 being a mixture distribution of uniform below LOD and normal above LOD, and the new Z 2 was used to simulate the outcome Y .In this scenario, we aim to investigate whether the LOD accommodation approaches hold when the distributions of exposures are different below and above LOD.In this example, we arbitarily assumed that the change point of distribution was exactly at LOD as a case study.In practice we may not know the change point unless there are external information.All the other settings are the same as in Scenario 1.
-Scenario 4. We assumed a null effect (i.e., β 2 = 0 ) of Z 2 for values below LOD and β 2 = 0.8 for val- ues above LOD.The other settings are the same as those in Scenario 1.This allows us to investigate whether the LOD accommodation approaches hold when the relationships between exposures and outcome are different below and above the LOD.Again, as a case study we arbitrarily picked the LOD as the changing point for simplicity, which may not happen in practice For each exposure and a given percent of values below LOD, we pre-determined the LOD values as the corresponding percentile from an independently simulated exposure dataset with sample size 20,000.With each simulated dataset, we first employed each of the five LOD accommodation approaches, then analyzed the resulting datasets with elastic net regression, WQS and BKMR under a unified formulation, that is, Y = h Z log + α T X + ǫ, where h Z log is the exposure- response function.Specifically, h Z log is β T Z log for elastic net, ψ(w T Z log ) for WQS with ψ being the total effect of a mixture, w as the vector of weights (or relative importance) and Z log as the pre-specified quantized Z log , and a general form h Z log for BKMR that allows non- linear relationship and interactions (see Appendix A).The R packages 'glmnet' [46], 'gWQS' [17] and 'bkmr' [47] with R version 4.2.1 (The R Foundation for Statistical Computing: http:// www.r-proje ct.org/) were used to implement these mixture methods.
In our implementations, all packages in R were applied as a default setting.Tuning parameters for elastic net were obtained from tenfold cross-validation.In WQS, we used quartiles of exposures after applying each approach for handling LOD with 200 bootstrap samples and 60% validation dataset.Five imputed datasets were generated for conventional and truncated MI approaches, and the final estimates of the MI and truncated MI were obtained using Rubin's rules [44].The R package 'bkmrhat' was used to combine the estimates of the MI and truncated MI in BKMR (https:// cran.r-project.org/ web/ packa ges/ bkmrh at/ index.html).We conducted 1000 simulation runs for each scenario.R code is available on GitHub at https:// github.com/ ml5977/ LOD_ accom modat ion.
The goal of our simulation study is to evaluate how different LOD accommodations influence the results of downstream mixture analysis.Note that since we simulated the data, all comparisons are made to estimates from the using the full datasets (i.e., not subject to LOD).We made this choice instead of comparing to the truth because some models are expected to exhibit biases even when all data are observed due to departure from the true underlying model, and certain model coefficients may have different interpretations.For example, elastic net regression explores a bias-variance trade-off, so we expect to see biases due to shrinkage [48].WQS is based on exposure quantiles, so all the parameters can be interpreted as the average effect when exposures increase by one quantile, whereas the parameter in the true underlying model represents the effect corresponding to a oneunit change.We also do not compare across the three mixture analysis methods, which is beyond the scope of the current study.
For elastic net regression and WQS, we reported the average bias and empirical standard error (SE) of the parameter estimations.For BKMR, using model assessment measures similar to those in Bobb et al. [11], we regressed the estimated exposure-response function h with each LOD accommodation approach on h from the full dataset and reported the average intercept, slope, R 2 , and standard error (SE) of h to assess the goodness of fit of the overall effects [11].An intercept close to 0 and slope close to 1 indicate no influence of the LOD accommodation approach on the downstream mixture analysis.We further reported posterior inclusion probabilities (PIPs) for each exposure.To be consistent with BKMR results in assessing overall effect, R 2 of regressing h from each LOD accommodation on h with the full dataset were also reported for elastic net regression and WQS.

NHANES data to explore the relationship between POPs and LTL
In addition to the simulation studies, we applied the above LOD accommodation approaches to the NHANES data collected between 2001 and 2002 as described in Gibson et al. [23] and Mitro et al. [32].We considered a subset of 1,003 participants who were over twenty years old, and provided blood samples and consented to DNA analysis, with sufficient stored samples to estimate telomere length, and without any missing values for individual exposures and covariates not related to LOD, as described in Gibson et al. [23].The Institutional Review Board of the National Center for Health Statistics approved the survey [49].
To be consistent with Gibson et al. [23], we restricted our analysis to 18 POPs with less than 40% of values below LOD, which include 11 polychlorinated biphenyls (PCBs), 3 dioxins, and 4 furans (Gibson et al. [23]).All samples were measured using high-resolution gas chromatography/isotope-dilution high-resolution mass spectrometry [50,51].LODs were typically ∼ 2ng/g , although they could be as high as 10.5ng/g [32], and 68.4% of sub- jects had at least one exposure below LOD.Using the data, Gibson et al. [23] and Mitro et al. [32] hypothesized that exposures to dioxins, furans, and PCBs were associated with longer LTL, which is the outcome of interest in this analysis.
We employed the above-mentioned approaches for accommodating exposures subject to LOD.All exposures were log-transformed due to their right-skewness.We adjusted for all the covariates as in Mitro et al. [32] and Gibson et al. [23], including age, age 2 , sex, race/ethnicity, educational attainment, BMI, serum cotinine, and blood cell count and distribution (white blood cell count, percent lymphocytes, percent monocytes, percent neutrophils, percent eosinophils and percent basophils).
Using the same data, Gibson et al. [23] handled exposure values below LOD through substituting them by LOD/ √ 2 .They found three potential drivers (PCB 126, PCB 118, and furan 2,3,4,7,8-pncdf ) selected by penalized regression methods, a positive overall effect of the POPs by WQS, a positive linear association with furan 2,3,4,7,8-pncdf, suggestive evidence of linear associations with PCBs 126 and 169, and a positive overall effect of the mixture but no interactions among exposures by BKMR.We re-analyzed the data with the same mixture methods after processing the values below LOD with five LOD accommodation approaches.
We recognized the need for sampling weights to account for the complex NHANES sampling scheme, in order to obtain results generalizable to the US population [49].However, our goal was to empirically compare the impact of different LOD accommodation approaches, rather than to provide estimates generalizable to the population.Thus, we simplified our analysis here by not including sampling weights for the NHANES cohort, so our results were consistent with those in Gibson et al. [23].We do recommend incorporating sampling weights into the analysis if an generalizable estimate is needed.We note some of the mixture analysis methods explored here, such as BKMR and WQS, require additional efforts to appropriately incorporate sampling weights, which is beyond the scope of this paper.

Simulation results: elastic net regression
Depending on the scenarios, the overall percent of subjects without any value below LOD in the simulated data was approximately 30% to 40%.Table 1 showed the bias of exposures Z 1 to Z 3 (group 1) and R 2 for each LOD accommodation approach with elastic net regression, while all other results for elastic net is in Table S1.
In Scenario 1 as a general case, when the exposures were moderately correlated, most approaches were unbiased except for the complete case analysis which also had higher SE, indicating inefficiency.In the high correlation setting, the biases in complete case analysis persisted, while imputing values below LOD by LOD/ √ 2 and con- ventional MI also showed biases for β 2 .The bias in MI decreased when truncated MI was used.Imputation by estimates from the AFT model (F-AFT) and truncated MI were empirically unbiased and efficient in both moderate and high correlation settings under Scenario 1.
When Z 2 was subject to 70% values below LOD and was completely ignored in the elastic net regression (Scenario 2A), all LOD accommodations performed poorly with low R 2 and large biases for exposures in the same group ( β 1 andβ 3 ) and covariate X 1 ( α 1 ).Note that exposures in other groups were relatively less impacted since the effect of Z 2 was potentially accounted for by those in the same group (i.e., Z 1 and Z 3 ).These biases decreased when correlations were higher, again presumably because the information in Z 2 was better captured by other exposures in the same group.These biases in β 1 and β 3 were further alleviated when an indicator variable of Z 2 (Scenario 2B), i.e., I(Z 2 > LOD ), was used.How- ever, β 2 now has a different interpretation (i.e., the differ- ence between values above LOD versus below the LOD), so we expected to see its large bias.Although Group 1 exposures were still biased in Scenario 2B for most of the LOD accommodation approaches, F-AFT and truncated MI generally performed well, especially in high correlation setting, followed by imputation by LOD/ √ 2 and conventional MI.
When different distributions below and above LOD were assumed for Z 2 (Scenario 3) or Z 2 had different effects below and above LOD (Scenario 4), all approaches for handling LOD, including model-based approaches such as truncated MI and F-AFT, performed poorly for β 2 because we lack any information to make inference about the values and relationship below LOD.Surprisingly, we observed a smaller bias of β 2 with conventional MI.However, the bias increased dramatically when the percent of values below LOD increased (results not shown).Furthermore, conventional MI was substantially biased in the intercept α 0 and was inefficient in β 2 com- pared to other LOD approaches with a lower R 2 .There- fore, truncated MI and F-AFT still performed relatively better than other approaches and using LOD/ √ 2 yielded slightly worse results but comparable.

Simulation results: WQS regression
We summarized the β 1 to β 3 and R 2 results in Table 2 and all remaining results for WQS in Table S2.We expected WQS to be less sensitive to values below LOD due to using quantized exposures.However, some LOD accommodations could disrupt the quantiles and result in large biases.For example, in Scenario 1, we found that F-AFT and truncated MI mostly maintained the exposures' quantiles and were empirically unbiased and efficient (Table 2).Complete case analysis showed relatively large biases, especially for overall effect estimate ( ψ ) in the setting of moderate correlation, due to the loss of all values below LOD and complete change of quantiles.Conventional MI also showed slightly larger biases compared to truncated MI because the imputed values could Table 1 Bias (SE) for exposures in Group 1 and R 2 with elastic net regression and each LOD accommodation approach compared to using full dataset Bias (SE) was reported for exposures in Group 1 ( β 1 , β 2 and β 3 ).All other results are provided in Table S1.All comparisons were made to the parameters with full datasets without LOD.R 2 was calculated by regression h from each LOD accommodation on h with the full dataset.In Scenario 2A, β 2 was not estimated because Z 2 was not included in the analysis Complete case -0.12 (0.26) -0.occasionally exceed the detection limit that can change quantile estimates.When LOD/ √ 2 was used, perfor- mance was unstable because the exposure's quantiles may not be maintained in the analysis of WQS if the percent of value below LOD is high (e.g., potential mis-assignment of quantiles).In evaluating the overall effects of the mixture with R 2 , complete case analysis underperformed across all LOD approaches while the others were similar.
When the percent of values below the LOD for Z 2 was increased to 70% and Z 2 was ignored in the WQS analy- sis (Scenario 2A), in the moderate correlation setting, the biases increased, especially for effects of exposures in the same group ( Z 1 and Z 3 ), total effect ψ , intercept and covariate X 1 .When an indicator variable I(Z 2 > LOD) was used as in Scenario 2B, the bias of total effects was slightly alleviated, but biases in weights of group 1 exposures, intercept and covariate X 1 persisted.All LOD accommodations performed similarly well in the high correlation setting, except complete case analysis was substantially biased in intercept and with lower R 2 .In the scenarios with different distributions (Scenario 3) or different effects (Scenario 4) below and above LOD for Z 2 , truncated MI and F-AFT maintained better perfor- mance in both parameter estimates and R 2 compared to Table 2 Bias (SE) for exposures in Group 1 and R 2 with WQS and each LOD accommodation approach compared to using full dataset Bias (SE) was reported for the total effect ( ψ ) and exposures in group 1 ( w 1 , w 2 and w 3 ).All other results are provided in Table S2.All comparisons were made to the parameters with full datasets without LOD.R 2 was calculated by regression h from each LOD accommodation on h with the full dataset.In Scenario 2A, w 2 was not estimated because Z 2 was not included in the analysis the other LOD accommodation approaches.Imputation by LOD/ √ 2 showed similar R 2 , but there was a large bias in estimating the total effect ψ.

Simulation results: BKMR
Table 3 showed the simulation results of BKMR under different scenarios.In Scenario 1, F-AFT performed the best among all the approaches, with intercept close to 0 and slope close to 1, indicating empirically unbiased results of h Z log .The F-AFT also led to high R 2 and lower SE.Truncated MI performed similarly to F-AFT but was slightly less efficient.Complete case analysis and imputation by LOD/ √ 2 underperformed, especially in the high correlation setting.In Scenario 2, none of the LOD accommodation approaches performed satisfactorily, despite the indicator variable (Scenario 2B) resulting in slightly better estimation than Scenario 2A.In Scenarios 3 and 4, F-AFT and truncated MI were the most unbiased and efficient in both correlation settings.In identifying important mixture components by PIPs, F-AFT and truncated MI performed similarly to using the full dataset, while complete case analysis showed discrepancies (Figure S2).The performance of imputation Table 3 Summary measures of estimated h(Z)with BKMR and each LOD accommodation approach compared to using full dataset Summary measures were obtained by regressing the estimated h of each LOD-accommodation approach on husing full datasets, and reported average intercept, slope, andR by LOD/ √ 2 in PIPs was comparable to those of F-AFT and truncated MI, despite this approach showing unstable results in some cases (e.g., high correlation settings).

NHANES data analysis results
When applying the elastic net regression to the mixture, F-AFT, truncated MI, and imputation by LOD/ √ 2 gen- erally resulted in similar findings (Fig. 1).Specifically, they all identified six important POPs: PCB99, PCB118, PCB126, PCB169, furan 2,3,4,7,8-pncdf, and furan 1,2,3,4,6,7-hxcdf with similar effects.Complete case analysis only identified PCB126 and PCB169 to be important, while conventional MI resulted in selecting many more exposures.We additionally conducted group lasso with the 18 POPs categorized into three groups: non-dioxinlike PCBs, non-ortho-PCBs, and mPFD, as described above.None of exposures in non-dioxin-like PCBs were selected except when using conventional MI, while nonortho PCBs (i.e., PCB126 and PCB169) were associated with non-zero coefficients in all LOD approaches (Figure

S3
).The magnitudes of the non-ortho PCB effects were larger with complete case analysis and conventional MI while the other three approaches yielded similar effects.For the mPFD exposures, again, F-AFT, truncated MI and imputaiton by LOD/ √ 2 estimated similar coeffi- cients and they all selected furan 2,3,4,7,8-pncdf as the most important exposures, followed by PCB118.Complete case analysis resulted in null effects for all mPFDs, and conventional MI showed mild effects for some nondioxin-like PCBs and opposite direction for some of the mPFDs.
Deciles in exposures were used in the analysis of WQS regression to be consistent with Gibson et al. [23].The total effect of 18 POPs ranged between 0.014 and 0.018 under various LOD-handling approaches, and they were statistically significant except for with complete case analysis, which is expected due to loss of efficiency with only 317 subjects included in the analysis (Table 4).Applying a priori cut-off weight of 1/18, we found 3 to 4 important POPs across these LOD accommodation approaches.Imputation by LOD/ √ 2 and truncated MI found 2,3,4,7,8-pncdf as the most important exposure, followed by PCB126 and 1,2,3,4,6,7,8-hxcdf.In addition to these three, the F-AFT approach also identified PCB194.
Using BKMR, we employed hierarchical variable selection with the three pre-defined groups, which provided importance scores for both the groups (i.e., group PIPs) and each exposure within a group (i.e., conditional PIPs).Truncated MI and imputation by LOD/ √ 2 both resulted in the non-ortho PCB group with the highest PIP among three groups, while mPFD group has the highest PIP with F-AFT, conventional MI and complete case analysis (Table S3).Within the mPFD exposures, furan 2,3,4,7,8pncdf contributed most to the model, followed by PCB 118 when imputation by LOD/ √ 2 , truncated MI and F-AFT approaches were used.PCB 169 and PCB 126 in the non-ortho PCB group had similar importance weights when we applied imputation by LOD/ √ 2 , trun- cated MI and F-AFT approaches.The individual effects of the POP exposures showed linear trends across LOD accommodation approaches (Fig. 2A), while the magnitudes of associations varied, especially for PCB169 and furan 2,3,4,7,8-pncdf which were selected as important exposures among the 18 POPs.The overall mixture effect was also close to a positive linear trend on the LTL outcome across LOD approaches, while the strength and efficiency varied (Fig. 2B).

Discussion
In this study we have compared how five popular approaches for handling exposures subject to LOD influence the results of mixture analysis.We did not mean to provide a guideline on how to handle values below LOD, rather to draw attention about how results can be misled by the various LOD accommodation approaches, and would like to advocate for careful examination of LOD accommodation prior to applying downstream mixtures analysis.
Through our extensive simulations, we generally favored using truncated MI and censored AFT models to impute values below LOD for the stability of downstream mixture analysis when the percent of the LOD was low to moderate (e.g., 30-50%).Compared to other approaches, truncated MI and censored AFT models generate imputed values based on the information from other exposures and covariates and guarantee that the imputed values are below LOD.Satisfactory results were also found with these two approaches when evaluating statistical uncertainties, such as mean squared error and coverage probability of the 95% confidence interval, in additional linear regression simulations (Table S4), as well as when incorporating grouping information in the analysis of group lasso and BKMR with hierarchical variable selection (Tables S5 and S6).Of course, these model-based approaches rely on modeling assumptions  and borrowing information from other exposures and baseline demographics.However, we argue since we do not get to observe any information below LOD, we need some assumptions, and the modeling assumptions made in these two approaches are relatively mild and reasonable in practice.
Complete case analysis and imputations by LOD/ √ 2 are frequently used in environmental health studies due to their easy implementation.However, we found that their performance can be quite unstable, especially in scenarios with high correlations or high percent of values below the LOD as commonly observed in environmental mixture studies.Richardson and Ciampi [38] also reported the bias in risk estimates when an arbitrary constant value such as LOD or LOD/2 was used to replace values below the LOD, and pointed out the magnitude of bias depends on the differences between the substitution value and true exposure distribution below LOD.
When the percent of values below LOD increased to 70% in our simulations, using an indicator variable of whether the values are above LOD performed better than excluding the exposure variable in the analysis.When other exposures were highly correlated with the exposure that had a high percent of values below LOD, its influence on the overall effect was limited because its information was well captured by other exposures.Based on our simulation studies with various percent of values below the LOD (results not shown), we recommend using the indicator variable approach when the percent of exposures below LOD is above 50%.For the NHANES data analysis, we restricted to exposures with less than 40% of values below LOD, in order to replicate the analysis in Gibson et al. [23].If we were to perform our own analysis, we will likely use 50% as a cutoff to include three additional POPs in the analysis.
We acknowledge that it is difficult to verify an assumed relationship or distribution between exposure subject to LOD and disease outcome for values below the LOD.To address this, we examined various LOD accommodation approaches assuming that the relationship has no impact below the LOD (Scenario 3) and the distribution is different for values below LOD (Scenario 4) as a case study.In our simulation study, none of the approaches for handling LOD in this study performed satisfactorily, which is similar to the results given by Ortega-Villa et al. [33] for a single exposure.In such cases, we recommend considering the binary indicator approach for exposures with suspected differential distribution or relationship with outcome [52], while truncated MI or F-AFT can still be used for all other exposures subject to LOD.Even though BKMR with imputation by LOD/ √ 2 , truncated MI and F-AFT performed satisfactorily in such scenarios due to its flexibility in allowing non-linear relationship, the implementation of the missing indicator approach could lead to further performance enhancement in BKMR.Yet, interpreting the estimated coefficient for the missing indicator within the exposure-response function of BKMR might prove challenging, especially when indicators are needed for multiple exposures.
We applied the LOD approaches to NHANES 2001-2002 where 18 selected POPs were subject to different proportions of values below the LOD.In our analysis, we did not include sampling weights because our goal was to understand the impact of different LOD accommodation approaches on downstream mixture analysis as a comparison with Gibson et al. [23], which did not incorporate sampling weights.To incorporate sampling weights, Zhang et al. [53] sampled one bootstrap sample with replacement from the NHANES data, with probabilities proportional to the sampling weights to test the results.We also implemented the same procedure.Although the mixture analysis results were different, we observed similar patterns across LOD accommodation approaches (results not shown).
In this study, we considered a two-stage approach as a practical implementation where we first performed the LOD accommodation to get a "complete" dataset, then conducted mixture analyses using this dataset.In the multiple imputation (MI) with or without truncation, we generated five imputed datasets, and combined the results of mixture analysis using the Rubin's rule [44], which takes imputation variability into account in the final results.However, single imputations by LOD/ √ 2 , and by estimates from the censored AFT model did not account for the uncertainty resulting from the imputation, which could lead to an overestimation of the precision.This can also be seen in the linear regression simulation results in Table S4, with somewhat worse coverage probabilities by F-AFT and LOD/ √ 2 .Neverthe- less, in our experience working with epidemiologists, this two-stage approach is highly preferred in practice due to its convenience.It requires handling the LOD only once and allows the resulting dataset to be used as the "true" dataset for multiple studies in the future.As mentioned above, Chen et al. [41] proposed a semiparametric multivariate AFT approach with multiple exposures to simultaneously model the exposures subject to LOD and the outcome, which accounts for uncertainty in the exposure assessment.This approach was applied to study the relationship between a panel of urinary trace metals and oxidative stress in pregnant women.The use of this powerful approach is limited by its computational complexity.Thus, it is of great interest to extend this approach to allow simultaneous modeling of the exposures subject to LOD with various mixture outcome models, and provide user-friendly software.Some analytical laboratories often provide the machine readings for specimens whose observed values is declared to be below the LOD, with the understanding that the specimen's level of analyte cannot reliably distinguished from zero; these readings may involve substantial measurements errors.Machine-read values have been often used in environmental mixture studies [54][55][56].However, we did not consider the machineread approach in our case study because it is difficult to justify the actual mechanism of the machine-raed approach given that each machine in each lab has its unique way of generating the reads, and the accuracy could vary dramatically.In our data analysis, NHANES 2001-2002 also did not provide machine-read values.
Here, we limited to three mixture analysis methods including elastic net regression, WQS, and BKMR which have been widely used in environmental mixture studies.We are aware of many other mixtures anslysis methods and performed simulations to understand the impact of LOD accommodations on these methods too.However, they were not included due to the length of the current manuscript.For example, Keil et al. [57] recently proposed a quantile-based g-computation method (q-gcomp) that builds up on WQS regression integrating its estimation procedure with a g-computation technique, which is widely used for causal inference [58].The q-gcomp method relaxes the unidirectionality and linearity assumptions of the WQS regression.Results were similar to those for WQS, which is likely due to their similar model structures and our simulated exposures were all in one direction (e.g., see Table S7 for q-gcomp results under Scenario 1).
Several methodological extensions are of interest for further exploration.First, in this study we assumed a linear combination of variables in MI and F-AFT for imputation.However, these approaches could allow nonlinearity and/or non-additivity for better recovering the true effects in the mixture setting.We also assumed all the effects were in the same direction with no interactions, which limits generalizability, and assumed a multivariate normal distribution for the exposures.Second, our study employed popular approaches for accommodating LOD before applying the mixture analysis methods to the revised data (i.e., a two-stage approach).Lastly, environmental mixture exposures are often repeatedly measured (i.e., longitudinal mixture exposures), which could allow more accurate modeling of the exposure trajectories.We leave a consideration of LOD adjustments that can appropriately incorporate longitudinal mixture exposures as a project for further development.

Conclusion
Quantifying the impact of mixtures of environmental exposures is becoming increasingly important for identifying disease risk factors and developing targeted public health interventions.Our case study delved into the issue of LOD in detail to understand how common approaches for handling LOD impact downstream mixture analysis.Our exploration provides insight into various LOD accommodation approaches in downstream mixture analyses, enhancing the quality and reliability of environmental health studies.

Appendix A. Detailed descriptions of mixture analysis methods
For each subject i(i = 1, . . ., n) , let Y i be a continuous outcome of interest.Let Z i and X i denote p -and q-vector of exposures and covariates, respectively.Note that X i includes 1 for the intercept term.Thus, we observe data {Y i , Z i , X i , i = 1, . . ., n} with sample size n .Using these notations, ordinary linear regression can be specified as Y i = β T Z i + α T X i + ǫ i , with ǫ i ∼ N 0, σ 2 .
Elastic net is a regularized regression method that incorporates the linear combination of L 1 and L 2 penalties of the lasso and ridge methods [12].The estimates can be obtained from argmin 2 �β� 2 2 + 2 �β� 1 , where tuning parameters 1 and 2 can be determined by cross validation (CV).Note that 2 = 0 and 1 yield ridge and lasso regressions, respectively.We used the R package 'glmnet' [46].
WQS regression can be specified as Y i = ψ w T Z i + α T X i + ǫ i , where Z is a pre-specified quantized vari- able of exposure Z .ψ represents the coefficient for the overall linear effect of the mixture, and w is the weight of each exposure.This method assumes that sum of all weights is 1 and each weight is between 0 and 1.Furthermore, this method assumes the same direction in all exposures (i.e., unidirectionality assumption).To conduct the WQS regression, we used 'gWQS' R package [17].
BKMR includes a completely nonparametric function of exposures as Y i = h Z 1i , . . ., Z pi + α T X i + ǫ i , where h(•) characterizes a high-dimensional expo- sure-response function that may incorporate nonlinearity and/or interaction among the mixture components.BKMR provides the posterior inclusion probabilities for each exposure, plotting the exposure-response function, and the cumulative (or overall) effects of the mixture.The R package 'bkmr' was used for the analysis [47].

Fig. 2
Fig. 2 Individual and overall relationships of 18 POPs with log-LTL from BKMR using NHANES 2001-2002 data.A Exposure-specific effect estimates of mixture members.B Overall effect of the mixture.Abbreviations: Imputation by LOD/ √ 2 (LOD/sqrt(2)); conventional multiple imputation (MI); truncated multiple imputation (Truncated MI); imputation by estimates using the AFT model (F-AFT)

Table 4
Total effect of 18 POPs and important exposures identified from WQS with each LOD accommodation approach, with NHANES 2001-2002Relative importance was reported for exposures with the weight over 1/18